Pfaffian pairing and backflow wave functions for electronic structure 

quantum Monte Carlo methods 
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We investigate pfaffian trial wave functions with singlet and triplet pair orbitals by quantum Monte 
Carlo methods. We present mathematical identities and the key algebraic properties necessary for 
efficient evaluation of pfaffians. Following upon our previous study,— we explore the possibilities 
of expanding the wave function in linear combinations of pfaffians. We observe that molecular 
systems require much larger expansions than atomic systems and linear combinations of a few 
pfaffians lead to rather small gains in correlation energy. We also test the wave function based on 
fully-antisymmetrized product of independent pair orbitals. Despite its seemingly large variational 
potential, we do not observe additional gains in correlation energy. We find that pfaffians lead 
to substantial improvements in fermion nodes when compared to Hartree-Fock wave functions and 
exhibit the minimal number of two nodal domains in agreement with recent results on fermion nodes 
topology. We analyze the nodal structure differences of Hartree-Fock, pfaffian and essentially exact 
large-scale configuration interaction wave functions. Finally, we combine the recently proposed form 
of backflow correlations^ with both determinantal and pfaffian based wave functions. 

PACS numbers: 02.70.Ss, 71.15.Nc, 71.15.-m 



INTRODUCTION 



One of the most promising many-body electronic struc- 
ture approaches is the quantum Monte Carlo (QMC) 
method, which employs stochastic techniques for solv- 
ing the stationary Schrodingcr equation and for evalua- 
tion of expectation values4£&£ QMC methodology has 
an important virtue that it enables us to test and em- 
ploy variety of many-body wave functions (WFs) with 
explicit electron-electron correlation. This opens a pos- 
sibility to explore wave functions, which are very difficult 
to use with traditional methods based on one-particle 
basis expansions and on orthogonality of one-particle or- 
bitals. These high accuracy wave functions enable us to 
understand the nature of many-body effects and also to 
decrease the QMC fixed-node errors which come from 
the necessity to circumvent the fermion sign problem. 8,9 
Fixed-node QMC has been very effective in providing 
high accuracy results for many real systems such as 
molecules, clusters and solids with hundreds of valence 
electrons. Typically, for cohesive and binding energies, 
band gaps, and other energy differences the agreement 
with experiments is within 1-3% The key challenge 
for successful application of fixed-node QMC is to de- 
velop methods, which can eliminate the fixed-node bias 
or at least make it smaller than experimental error bars 
for the given quantity. This is a difficult task, once we re- 
alize that the fermion nodes (the subset of position space 
where the wave function vanishes) arc complicated high- 
dimensional manifolds determined by the many-body ef- 
fects. So far, improvement in the accuracy of trial wave 
functions has proved to be one realistic approach to find- 



ing better approximations for the nodes. This approach 
has an additional benefit in forcing us to think about the 
relevant correlation effects and their compact and com- 
putationally efficient description. 

The commonly used QMC trial wave functions have 
the Slater-Jastrow form, which can be written as = 
exp[L r corr ] where ^a is the antisymmetric part while 
Ucorr describes the electron-electron and higher-order 
correlations. The antisymmetric component is typically 
one or a linear combination of several Slater determi- 
nants of one-particle orbitals, such as Configuration In- 
teraction (CI) expansion. 11 To overcome the limit of one- 
particle orbitals, the two-particle or pair orbital has been 
suggested. In condensed systems, one such example is 
the Bardeen-Cooper-Schrieffer (BCS) wave function, 12 
which is an antisymmetrized product of identical singlet 
pairs. 

The pair orbital is also referred to as geminal and the 
resulting wave function as the antisymmetrized product 
of geminals (APG) i 13 ' 14 In its full variational limit, say, 
for a system with M pairs of fermions, such a wave func- 
tion is an antisymmetrized product of M distinct pair 
orbitals. However this freedom leads to computationally 
very demanding form, since the number of terms grows 
exponentially with the number of pairs. In this work, we 
have actually tested the APG wave functions on a few- 
particle systems as will be reported later. On the other 
hand, if all the geminals in the product are identical, then 
this special case of APG is called the antisymmetrized 
geminal power (AGP). — It can be shown that for a sin- 
glet type geminal with a fixed number of particles the 
AGP and BCS wave functions are identical— and can be 
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expressed in the form of single determinant 

The BCS wave function has been recently used to 
calculate several atoms and molecules as well as su- 
perfluid Fermi gases i 18 i 19 ' 20 The results show promising 
gains when compared to the single-determinant Hartree- 
Fock (HF) wave functions, nevertheless, in partially spin- 
polarized systems the improvements are less pronounced 
due to the lack of pair correlations in the spin-polarized 
subspace . 18 ' 19 The spin-polarized (triplet) pairing wave 
functions naturally lead to pfafhans (instead of determi- 
nants). In this respect, pfafhan have been mentioned 
a few times and applied to model system o 21 i 22 ' 23 in the 
past. 

In this paper, we follow upon our previous letter^ in 
which we have proposed the description of electron sys- 
tems by pfaffian wave functions with variational freedom 
beyond HF and BCS wave functions. The pfaffian form 
proved to be the right algebraic form, which can accom- 
modate pair orbitals with singlet and triplet pair chan- 
nels, together with unpaired one-particle orbitals, into a 
single compact wave function. Here we present a set of 
key mathematical identities and formulas for pfafhans, 
some of them derived for the first time. We investigate 
generalizations to linear combinations of pfaffians and to 
antisymmetrized independent singlet pairs and compare 
the results from the point of view of recovered energies 
and compactness of the wave functions. We analyze the 
fcrmion nodes for some of the wave functions and point 
out the topological differences between HF, pfaffian and 
essentially exact wave functions for a given test exam- 
ple. Finally, we explore the possibility of further im- 
provements of nodal structure of pfaffians by using the 
recently proposed form of backflow correlations^ 



II. ALGEBRA OF PFAFFIANS 

A. Definitions 

First introduced by Arthur Cayley in 1852)21 the 
pfaffian is named after German mathematician Johann 
Friedrich Pfaff. Given a 2n x 2n skew-symmetric matrix 
A 



the pfaffian of A is defined as antisymmetrized 



product 



pf [A] = A[ai^a 3A . . . a 2n -i.2n] 



X " sgn(a) 

a ii ,ji a i 



(1) 



where the sum runs over all possible (2n — 1)!! pair parti- 
tions a = {(ix,jx), (i 2 , 32), (injin))} of {1, 2, .. . ,2n} 
with ik < jk- The sign of permutation associated with 
the partition a is denoted as sgn(a). The pfaffian for a 
matrix of odd order equals to zero. The following exam- 



ple gives pfaffian of a A(A x 4) skew-symmetric matrix: 
pf 



an 

-an 

-ai3 — 0,23 

-<2i4 — a 2 4 



013 ai4 

«23 a 24 

a 34 

-»34 . 



»12«34 — ^13^24 + ai4fl23- 

(2) 



It can be also evaluated recursively as 



j = 2 Ql.j 

2n 

3=2 



(3) 



where ttxj is partition with ik,jk 7^ l)j an d P c (ax,j) is 
defined as pfaffian cofactor of ai.j. The cofactor for an 
element a,j,k is given by a formula 



P c (a jtk ) = (-iy+ k+1 pf[A(j,k;j,k)], 



(4) 



where the matrix A(j, k; j, k) has the rank 2(n— 1) x 2(n— 
1) and is obtained from A by eliminating rows j and k 
and columns j and k. 



B. Calculation of a pfaffian 

There exist several identities involving pfaffians and 
determinants. For any 2n x 2n skew-symmetric matrix 
A and arbitrary matrices B(2n x 2n) and M(n x n) we 
have the following relations: 



pf 



pf 



P fL4 T ] 

P fL4] 2 

Ax 
A 2 

vi[BAB T 

M 
-M T 



(-l)>f[A] (5a) 

det[A] (5b) 

pf[Ai]pf[A 2 ] (5c) 

det[B]pf[A] (5d) 

(-l) n ^det[M] (5e) 



Proofs: 

(|5ap Each permutation contains product of n pairs re- 
sulting in an overall (—1)™ factor. 

(|5b[) This is a well-known Cayley's relationship be- 
tween the pfaffian and the determinant of a skew- 
symmetric matrix. Since it has been proved many times 
before in variety of wayS ) 25 ' 26 i 27 we do not give this proof 
here. Using this relation we rather prove a more gen- 
eral version of Cayley's identity 2 "^ in the Appendix (|A"|) . 
which we were not able to find anywhere else except in 
the original Cayley's paper.— 

(|5cl) Use the expansion by pfaffian cofactors. 

(|5d[) By squaring (4d), using Eq. (|5b[) , and taking the 
square root one finds pf[BAB T ] = ±det[B]pf [A], Sub- 
stituting the identity matrix I for B one finds + to be 
the correct sign. 
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(l5el) Assume 



B 



M 
I 



and A 



I 

-I 



in Eq. (|5d|). The overall sign is given by value of pfL4]. 

The identities listed above imply several important 
properties. First, Eqs. (|5d|) and (|5e|) show that every 
determinant can be written as a pfaffian, but on the con- 
trary, only the absolute value of pfaffian can be given by 
determinant [Eq. (|5b[) ] . The pfaffian is therefore a gener- 
alized form of the determinant. Second, by substituting 
suitable matrices 28 for M in Eq. (|5d|) one can verify the 
following three properties of pfaffiansj 2 ^ similar to the 
well-known properties of determinant: 

(a) multiplication of a row and a column by a con- 
stant is equivalent to multiplication of pfaffian by 
the same constant, 

(b) simultaneous interchange of two different rows and 
corresponding columns changes the sign of pfaffian, 

(c) a multiple of a row and corresponding column 
added to to another row and corresponding column 
does not change the value of pfaffian. 



Any real skew-symmetric matrix can be brought to block- 
diagonal form by an orthogonal transformation. Recur- 
sive evaluation [Eq. ([3])] then implies that the pfaffian of 
block-diagonal matrix is directly given by 



pf 



Ai 
-Xi 



A 2 
-A 2 



A„ 

-An 0. 



AiA 2 ...A„. (6) 



Therefore by employing a simple Gaussian elimination 
technique with row pivoting we can transform any skew- 
symmetric matrix into block-diagonal form and obtain 
its pfaffian value in 0(n 3 ) time. 

However, in QMC applications one often needs to eval- 
uate the wave function after a single electron update. 
Since Cayley 2 ^ showed [for proof see App. (TAJ] that 
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— CJ3,n • 


. 



= pf 






ai2 


£113 


• ■ ai,n 







&12 


bis • 


• bi t „ 


— ai2 





CI23 


• ■ £12, n 




— b V 2 





ffl23 


■ ■ a,2,n 


-ai3 


-C!23 





• ■ £t3, n 


Pf 
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we can relate the pfaffian of original matrix pf [A] to the 
pfaffian of a matrix with updated first row and column 
pf[-B] using the inverse matrix A -1 in only 0(n) opera- 
tions by 



pf[B] = 



det[A] EjhjAj, 1 
pf[A] 



= P f[A]J2b lj Aj 1 1 . (8) 



The second part of Eq. ([5]) was obtained by taking ad- 
vantage of the identity in Eq. (|5b|) . Slightly more com- 
plicated relation between pi[A] and pf[-B] can be de- 
rived if one considers simultaneous change of two separate 
rows/columns, which represents a two electron update of 
the wave function. 



fian: 



1 gpfjA] _ 1 
pf[A] dci ~ 2 1 



(9) 



and 

1 <9 2 pfL4] _1 
pfp] dc.dcj ~2 



tr 



1 

-tr 

4 



_! C9 2 A ' 




1 

-tr 
2 


A- 1 


dcidcj 




dc l _ 


tr 







-A' 



where A 1 is again the inverse of A. 



(10) 



C. Gradient and Hessian of pfaffian 

If the elements of matrix A depend on some variational 
parameters {c^}, one can derive the following useful ex- 
pressions (see Sec. IIV|) for gradient and Hessian of pfaf- 



III. 



PAIRING WAVE FUNCTIONS 



In order to contrast the properties of pair wave func- 
tions with the wave functions build from one-particle or- 
bitals we will first recall the well-known fact from the 
Hartree-Fock theory. The simplest antisymmetric wave 
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function for N electrons constructed from one-particle 
orbitals is the Slater determinant 



^hf = det[ip k (n, s^} = det[<£k(i)]; i,k = l, 



,N, 
(11) 



where tilde means that the one-particle states depend on 
both space and spin variables. Clearly, for N electrons, 
this requires N linearly independent spin-orbitals which 
form an orthogonal set in canonical HF formulation. 

Let us now consider the generalization of the one- 
particle orbital to a two-particle (or pair) orbital <j>(i,j), 
where tilde again denotes dependence on both spatial and 
spin variables. The simplest antisymmetric wave func- 
tion for 2N electrons constructed from the pair orbital is 
a pfaffian 



* = ^(1,2),^(3,4)...^(2JV 



l,2AT)]=pffe)]- 
(12) 



The antisymmetry is guaranteed by the definition l|T]). 
since the signs of pair partitions alternate depending on 
the parity of the corresponding permutation. The im- 
portant difference from Slater determinant is that in the 
simplest case only one pair orbital is necessary. (This can 
be generalized, of course, as will be shown later.) If we 
further restrict our description to systems with collinear 
spins, the pair orbital <p{Ti, rj, Sj) for two electrons in 
positions and rj and with spins projections Sj and Sj 
and can be expressed as 

fo i ,s i ;T j ,s j )=<l>(i,j){8 i 8 j \[\ T|> - I 4T}]/\/2 (13) 

+ x 1 Hhj)(s i s j \ TT) 

+ X U (iJ)(siSj\[\ U) + UT>]/V2 
+ X U (h.i)(stS 3 \ ||). 

Here <fi(i,j) = 4>(ri,Tj) is even, while x > X an( l X 
are odd functions of spatial coordinates. In the rest of 
this section we will discuss special cases of wave function 

an). 



A. Singlet pairing wave function 

Let us consider the first i — 1, 2, TV electrons to be 
spin- up and the rest j = N + 1, ...,2N electrons to be 
spin-down and allow only </>(r;,rj) in </>(rj, Si\rj, Sj) to 
be non-zero. Using the pfaffian identity [Eq. (|5e))] we can 
write the wave function for N singlet pairs, also known 
as the BCS wave function, in the following form 



BCS 



Pf 



* n 

-*TIT 



det[* u ] 



(14) 



which is simply a determinant of the NxN matrix = 
[(j>(i,j)] as was shown previously J^- 

It is straightforward to show that the BCS wave func- 
tion contains the restricted HF wave function as a special 



case. Let us define the Slater matrix C — [<Pi(j)], where 
{ifii} is a set of HF occupied orbitals. Then we can write 

V HF = det[C]det[C] = det[CC T ] = det[*# F ], (15) 
where 



a>TJ- \ 

HF> 



<t>HF(i,j) 



N 

E 



<Pk(i)<Pk(j)- (16) 



On the other hand, we can think of the BCS wave func- 
tion as a natural generalization of the HF one. To do so 
we write the singlet pair orbital as 



>N 



) =^2S k ,i<pk(i)<Pi(j) = (p(i)S(p(j), 



(17) 



k.l 



where the sum runs over all (occupied and virtual) single- 
particle orbitals and S is some symmetric matrix. There- 
fore, we can define one-particle orbitals which diagonal- 
ize this matrix and call them natural orbitals of a singlet 
pair. 

The BCS wave function is efficient for describing sys- 
tems with single-band correlations such as Cooper pairs 
in conventional BCS superconductors where pairs form 
from one-particle states close to the Fermi level. 



B. Triplet pairing wave function 

Let us assume, in our system of 2N electrons, that the 
first Mi are spin-up and remaining — 2N — Mi are 
spin-down. Further we restrict Mi and Mi to be even 
numbers. Then by allowing only \ and x (hj) 
in (fl3|) to be non-zero, we obtain from (fT2|) by the use of 
Eq. (09 



* 7 



pf 



£TT 
£ U 



pf[£ TT ]pfK 



111 



(18) 



where we have introduced Mi x Mi(M2 x M2) matrices 
£Tt(U) = r^Tt(U) (j 5 j)l _ To our knowledge, this result 
was never explicitly stated and only the weaker statement 
that the square of wave function simplifies to a product 
of determinants has been given^i 

The connection to a restricted HF wave function for 
the above state can be again established as follows. In 
accord with what we defined above, det[C^^] are spin- 
up(-down) Slater determinants of some HF orbitals {tfi}. 
Then, by taking advantage of Eq. (|5e|) we can write 



$>hf = det[C T ]det[C i ] 



(19) 



pf[c T AiCH P f[cM 2 c+ 

pf[Ai]pf[A 2 ] 



given Ai and A2 are some skew-symmetric non-singular 
matrices. In the simplest case, when A\ and Ai have 
block-diagonal form with all values \ = 1, one gets 



^hf = P f [€£ F ]pf ie H i F ] 



(20) 
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The pair orbitals can be then expressed as 



(21) 



Mi(M 2 )/2 

E 

k=l 



[<P2k-l(i)<P2k(j) ~ <P2k-l(j)<P2k(i)]- 



Similarly to the singlet pairing case, one can also think 
of triplet pairing as a natural generalization of the HF 
wave function. To do so we write the triplet pair orbitals 

as 

>Mi(M 3 ) 

x(ij) mi) = E 4!/ u V*(i)wC7) 

k.l 



(22) 



where again the sum runs over all (occupied and vir- 
tual) single-particle orbitals and A^TCii) are som e skew- 
symmetric matrices. Therefore, we can define one- 
particle orbitals which block-diagonalize these matri- 
ces and call them natural orbitals of a triplet spin-up- 
up (down- down) pair. 



pfaffian pairing wave functions. As we mentioned ear- 
lier, the trial variational wave function is a product of 
an antisymmetric part \& a times the Jastrow correlation 
factor 

tf T (R) = * a(R) exp[U corr ({r tJ }, {r u }, {r u })}, (24) 

where U corr depends on electron-electron, electron-ion 
and, possibly, on electron-electron-ion combinations of 
distancesi 7 ' 30 i 31 Previously, we have reported^ the results 
with antisymmetric part being equal to ^a = ^hf, 
^ a = ^bcs and ^a — ^stu- We extend this work 
to different linear combinations of pfaffians. The pair 
orbitals were expanded in products of a one-particle or- 
bital basis according to Eqs. (fT71 |2"2"|) . The expansions 
include both occupied and virtual one-particle orbitals 
from either Hartree-Fock or CI correlated calculations. 11 
The pair orbital expansion coefficients were then opti- 
mized in VMC by minimizations of energy, variance or 
a combination of energy and variance^ The optimiza- 
tion procedure requires the calculation of gradient and 
the Hessian of the wave function according to Eqs. ([HJ 
ITUl) . We used pseudopotentials^ to eliminate the atomic 



C. Generalized pairing wave function 



A. Single- and multi-pfaffian calculations 



Let us now consider a partially spin-polarized system 
with unpaired electrons. In order to introduce both types 
of pairing we allow x (h j)i X (i>j) an d 4>(i,j) m tHI 
to be non-zero. However, we omit the X (hj) term. 
Then our usual ordered choice of electron labels, with 
all spin-up electrons first and remaining electrons spin- 
down, enables us to directly write from (fP2f the singlet- 
triplet-unpaired (STU) orbital pfaffian wave function 1 



STU 



Pf 



(23) 



where the bold symbols are block matrices and vectors of 
corresponding orbitals as defined in sections IIII Al IIIIBI 
and T denotes transposition. Let us note that for a spin- 
restricted STU wave function the pair and one-particle 
orbitals of spin-up and -down channels would be identi- 
cal. 

The above pfaffian form can accommodate both singlet 
and triplet pairs as well as one-particle unpaired orbitals 
into a single, compact wave function. The correspon- 
dence of STU pfaffian wave function to HF wave function 
can be established in a similar way to the pure singlet and 
triplet pairing cases. 



The results of a single-pfaffian wave function calcula- 
tions when applied to the first row atoms and dimers 
were reported in our previous paper. 1 However, for com- 
pleteness, we summarize all total and correlation energies 
in Table IIIII Appendix [Bj A systematic high percentage 
of recovered correlation energy on the level of 94-97% 
in DMC method with generally low triplet contributions 
was observed. 

Further, the tests of multi-pfaffian (MPF) wave func- 
tion of the form 



y M PF=pf\Xl ,Xl, 01, Vl] +P%2 ,X2 » 4>2,<P2]+ ... 

(25) 

for atomic systems were also discussed in our previous 
study. 1 A small number of pfaffians was shown to recover 
another significant fraction of the missing correlation en- 
ergy comparable to much more extensive configuration 
interaction expansions in determinants. 

In this work, we extend the application of MPF wave 
functions to the diatomic cases (Table H]). However, 
only very limited gain over single STU pfaffian WF was 
achieved for MPF wave functions with few pfaffians. We 
therefore conclude that for obtaining significantly larger 
gains in correlation energy the molecular wave functions 
require much larger expansions. 



IV. PAIRING WAVE FUNCTION RESULTS 



APG wave function 



We perform the variational and fixed-node diffusion 
Monte Carlo (VMC and DMC) calculations^ with the 



We have also tested the fully antisymmetrized product 
of geminals (or singlet independent pairs) wave function, 
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TABLE I: Total energies for C and N dimers with amounts of correlation energy recovered in VMC and DMC methods with 
wave functions as discussed in the text. The corresponding number of pfaffians/determinants n for each wave function is also 
shown. See caption of Table 11111 for a more detailed description. 



Method 


WF 


n 


C 2 


E CO rr[%] 


n 


N 2 


E CO rr [%] 


VMC 


MPF 


5 


-11.0187(2) 


87.8(1) 


5 


-19.8357(3) 


89.2(1) 




APG 


4! 


-11.0205(4) 


88.3(1) 


5! 


-19.8350(3) 


89.0(1) 




CI" 


148 


-11.0427(1) 


93.7(1) 


143 


-19.8463(9) 


91.6(2) 


DMC 


MPF 


5 


-11.0437(4) 


94.0(1) 


5 


-19.8623(5) 


95.3(1) 




APG 


4! 


-11.0435(7) 


94.0(2) 


5! 


-19.8611(3) 


95.0(1) 




CF 


148 


-11.0573(2) 


97.3(1) 


143 


-19.875(2) 


98.3(5) 



a The determinantal weights were re-optimized in the VMC 
method. 



which introduces one pair orbital per each electron pair. 
For system of 27V fermions in singlet state, the APG wave 
function can be written as 

*4pg - ^i(l, 2), &(3, 4), . . . , $ N (2N - 1, 27V)] (26) 

= Pf ih (1, 2), 4 (3,4),..., 4, (27V - 1, 2N)], 
p 

where the last equation corresponds to the sum over all 
TV! possible permutations of TV different pair orbitals <j>i 
for each pfaffian. Recently, the APG wave function was 
used also by Rassolov^ 3 . in the form of an strongly or- 
thogonal geminals. Our results for C and TV dimers using 
APG wave functions in the VMC and DMC methods are 
given in Table HI 

Consideration of independent pairs results in an expo- 
nential increase of number of pfaffians. However, cap- 
tured correlation energy is on the level of small MPF ex- 
pansion and significantly less than CI with re-optimized 
weights using the same one-particle orbitals. This sug- 
gests that to achieve more correlation energy in larger 
systems, we have to go beyond double pairing. 

C. Nodal properties 

The fermion node is defined by an implicit equa- 
tion ^(R) — and for TV electrons it is a 
(37V — l)-dimensional hypersurface. With exception of 
few exact cases, the nodes of trial/variational wave func- 
tions introduce bias into fixed-node DMC energies. Re- 
cently a number of authors have reported improvement 
in nodal structure of trial wave functions ^^ i 18 ' 19 i 34 ' 35 

The effect of pairing correlations on nodes can be high- 
lighted by direct comparison. Fig. [T] shows the exam- 
ple of nodal structure of oxygen atom. Here we com- 
pare the nodal surfaces of HF (no pairing), MPF pfaffian 
(STU pairing) and a high accuracy CI wave function with 
more than 3000 determinants, which gives essentially ex- 
act fermion nodes [i.e., 99.8(3)% of correlation energy in 
fixed-node DMC]. 

It is clear that the changes in the nodal surfaces are 
significant, the most important one being the elimina- 
tion of artificial four nodal cells resulting from the in- 



dependence of spin-up and -down channels in HF. The 
pfaffian smooths-out the crossings and fuses the com- 
partments of the same sign into the single ones. These 
topology changes therefore lead to the minimal number of 
two nodal cells, an effect observed in correlated context 
previously i 36 i 37 ' 38 i 39 However, the nodes of the pfaffian 
wave functions could be further improved if the scheme 
for direct optimization of nodes of trial wave functions 
were usedi 20 ' 40 Additional result from our work is that 
despite such a substantial change in the nodal structure 
the amount of missing correlation energy is still non- 
negligible. 



V. BACKFLOW CORRELATED WAVE 
FUNCTIONS 

Another route to improvement of the trial wave func- 
tion and its nodal structure is through the introduction of 
backflow correlations^iii^ii 3 .^^^^^^ Given the 
form of our trial wave function, Eq. (|2~4")) . its nodal struc- 
ture is completely defined by the nodes of the anti- 
symmetric part <J/yi(R). The backflow correlations are 
then introduced by replacing ^(R) by ^(X), where 
X = (xi,X2, . . .) are some quasi-coordinates dependent 
on all electron positions R, such that overall antisym- 
metry is preserved. Consequently, if X is made depen- 
dent on some variational parameters, we can decrease the 
fixed- node errors by further optimizing ^^(X). 

The implementation of the backflow correlations into 
Slater determinant and pfaffian wave functions closely 
follows the approach of Kwon et al£& and Rios et al& 
The quasi-coordinate of i-th electron at position is 
given as 

Xj = Vi + &(R) 

= r 4 + $f»(R) + £f (R) + £f"(R), (27) 

where £j is the i-th electron's backflow displacement 
divided to the contributions from one-body (electron- 
nucleus), two-body (electron-electron) and three-body 
(electron-electron-nucleus) terms. They can be further 
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FIG. 1: (color online) A three-dimensional cut through the fermion node hypersurface of oxygen atom obtained by scanning 
the wave function with a spin-up and -down (singlet) pair of electrons at equal positions, while keeping the rest of electrons 
at a given VMC snapshot positions (small green spheres). Nucleus is depicted in the center of the cube by the blue sphere. 
The three colors (from left to right) show nodes of: Hartree-Fock (red/dark gray); multi-pfaffian nodes (orange/medium gray); 
and the nodes of the CI wave function (yellow/light gray) in two different views (upper and lower rows). The CI nodal surface 
is very close to the exact one (see text). The HF node clearly divides the space into four nodal cells while pfafnan and CI 
wave functions partitioning leads to the minimal number of two nodal cells. The changes in the nodal topology occur on the 
appreciable spatial scale of the order of 1 a.u. 



expressed as 



CT(R) = E^) r ^ 



(28) 



where = r,; 



Tn = T"i — Tj and we sum over all 
nuclei / and electrons j. The rj and 0\ with 9 2 terms 
are similar to one, two and three-body Jastrow terms 
present in U corr of trial wave function, Eq. (|24[) . and are 



further expanded as 

K r ) = ^2c k a k (r), 

k 
k 



kirn. 



VI. BACKFLOW WAVE FUNCTION RESULTS 

In this section, we present VMC and DMC results ob- 
tained with above implementation of backflow correla- 
tions for determinant and pfaffian wave functions. The 
Jastrow factors and pseudopotentials are identical to ones 
used in Sec. IIVI The distance-dependent basis func- 
tions {a} and {b} used in the Eq. (f!?9"|) are chosen ei- 
ther as Gaussians centered on the nucleus or polynomial 
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TABLE II: VMC and fixed-node DMC energies per electron 
and variances of local energies for various trial wave functions 
(S, Slater; SJ, Slater-Jastrow; SJBF, backflow correlated SJ) 
for 3D unpolarized HEG of 54 electrons. 





Method 


WF 


S/iV[a.u./electron] 


a 2 [a.u. 2 ] 


1.0 


HF 


S 


0.56925(2) 


19.3(1) 




VMC 


SJ 


0.53360(4) 


1.26(4) 






SJBF 


0.53139(4) 


0.81(4) 




DMC 


SJ 


0.53087(4) 


- 






SJBF 


0.52990(4) 


- 


5.0 


HF 


S 


-0.056297(7) 


0.776(4) 




VMC 


SJ 


-0.075941(6) 


0.0636(1) 






SJBF 


-0.078087(4) 


0.0236(1) 




DMC 


SJ 


-0.07862(1) 








SJBF 


-0.07886(1) 




20.0 


HF 


S 


-0.022051(2) 


0.0482(1) 




VMC 


SJ 


-0.031692(2) 


0.000711(4) 






SJBF 


-0.031862(1) 


0.000437(1) 




DMC 


SJ 


-0.031948(2) 








SJBF 


-0.032007(2) 





Pade functions 3 -^ to preserve the electron-electron and 
electron-nucleus (when used with pseudopotentials) cusp 
conditions ! 50 ! 51 The sets of variational parameters {c}, 
{d} and {g} are minimized in the similar fashion as in 
Sec. IIVI with respect to energy or mixture of energy and 
variance^ - In addition, all electron-electron coefficients 

1(2) 

{{dk} and {g^im} with fixed k and I) are allowed to be 
different for spin- like and for spin- unlike electron pairs. 



A. Homogeneous electron gas 

We benchmark our implementation of the backflow 
correlations on the homogeneous electron gas (HEG). 
The HEG system of 54 unpolarized electrons in the sim- 
ple cubic simulation cell with periodic boundary condi- 
tions was studied before^ ' 48 ! 49 We use the backflow dis- 
placement, Eq. P?)) , with only £f e being non-zero and 
let ry(r) in Eq. (|28p to be different for spin-like and for 
spin-unlike electron pairs. These functions are further ex- 
panded in the basis of polynomial Pade functions 3 -^ with 
cutoff equal to half of the simulation cell. 

We compare our results for the following three densities 
of r s = 1, 5 and 20 (see Table HI)) . First, it is clear that the 
HF and Slater-Jastrow (SJ) fixed-node DMC energies 
are in good agreement with previous results ^ 48 ' 49 Sec- 
ond, due to the omission of the three-body correlations 
from Jastrow factor and also from backflow displacement, 
it is expected that we obtain higher VMC energies and 
variances for SJ and backflow displaced SJ (SJBF) trial 
wave functions. Nevertheless, our fixed-node DMC ener- 
gies for SJBF trial wave functions closely match the re- 
sults of Kwon et al.^ and only slightly deviate at higher 
densities from results od Rios et al£ 



B. Carbon atom and dimer 

The backflow correlations in single-determinant 
Slater-Jastrow trial wave functions were recently applied 
also to inhomogeneous systems They were demon- 
strated to capture additional few percent of the correla- 
tion energy but being somewhat shy of the goal of more 
than 99%, with the only exception of Li atom. It was also 
suggested that the backflow by itself is unlikely to change 
the number of nodal cells. These observations let to fur- 
ther studies of backflow combined with the wave func- 
tions that have the minimal number of nodal cells — an 
important topological property associated with ground 
state wave functions. 36 ' 37 ! 38 ! 39 One of the successes of 
this scheme are the very recent results of Brown et alM- 
obtained from backflow correlated Cl-Jastrow wave func- 
tions applied to first row all-electron atoms. In this study, 
we further test the limits of the backflow correlations to 
decrease the fixed node errors of the Cl-Jastrow wave 
functions and, for the first time, also include backflow 
into the pfaffian-Jastrow pairing wave functions. Below 
is a brief discussion of our implementation and results for 
carbon atom and dimer systems. 

In the inhomogeneous backflow, each electron's co- 
ordinate is correlated by the displacement as given by 
Eq. (P?)) . while the functions rj and 6*1,2 are allowed to be 
spin-dependent. We use up to 11 Gaussian basis func- 
tions to fit the fi and rj functions, while the three-body 
functions 6*1,2 are limited to a product of 4 x 4 x 4 Gaus- 
sians. The main results are plotted in Figs. [2] and [3] and 
detailed numerical results are summarized in Tables IIVI 
and [V] of Appendix [Bj The backflow correlations are able 
to capture additional few percent of the correlation en- 
ergy for both Slater-Jastrow and pfaffian-Jastrow wave 
functions. Another important feature of backflow is 20- 
30% decrase in variances of local energy with respect to 
the wave functions without backflow correlations. We 
find that for the fully optimized backflow, the spin-unlike 
electron-electron functions are almost order of magnitude 
larger than spin-like ones as well as electron-nucleus func- 
tions. The gains are systematic with increasing number 
of parameters, however we do not find the three-body 
terms as important as reported in previous study 1^ This 
difference can be attributed to two main reasons — we use 
a different basis to expand the three-body functions 6*1.2 
and we also eliminate atomic cores by pseudopotentials. 
It is plausible that for systems with core electrons the 
three-body correlations are more important due to the 
strong variations of orbitals close to the nucleus. 

Finally, let us discuss the difference between the two 
systems with respect to missing correlation energy. For 
the C atom, we have shown previously^ that less than 
100 determinants gives more than 99% of correlation 
energy (E corr ). The C dimer 's fixed node errors are 
more pronounced, since the 148 determinants with re- 
optimized weights give only 97.5(1)% in a close agree- 
ment with recent calculations by Umrigar et al.£^ Em- 
ploying backflow correlations for our 148 determinant CI- 
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FIG. 2: (color online) Upper figure: Percentages of corre- 
lation energy from VMC and DMC methods versus a num- 
ber of backflow parameters with backflow correlated Slater- 
Jastrow (SJBF), pfaffian-Jastrow (PFBF) and Cl-Jastrow 
(CIBF) trial wave functions for C atom (2B, electron- nucleus 
and electron-electron terms; 3B, all electron-electron-nucleus 
terms; 23B for all terms together). Lines connecting the 
points serve only as a guide to the eye. Lower figure: Variance 
of the local energy versus a number of backflow parameters. 

Jastrow wave function gives no apparent gain in E corr ex- 
cept for decrease in the variance of local energy. The im- 
provement for the pfaffian-Jastrow wave function is also 
very modest (less than 1%). Our results suggest that to 
reach beyond 99% of correlation one still needs compli- 
cated multi-reference wave functions, even after including 
quite general forms of the backflow correlations. 

VII. CONCLUSIONS 

To summarize, we have proposed pfaffians with singlet 
pair, triplet pair and unpaired orbitals as variationally 
rich and compact wave functions. They offer significant 
and systematic improvements over commonly used Slater 
determinant-based wave functions. We have included a 
set of key mathematical identities with proofs, which are 
needed for the evaluation and update of the pfaffians. We 
have also shown connections of HF and BCS wave func- 
tions to more general pfaffian wave function. Further, 
we have explored multi-pfaffian wave functions, which 
enabled us to capture additional correlation. While for 
atomic systems the results are comparable to large-scale 
CI wave functions, 1 equivalent accuracy in molecular sys- 
tems most probably require much larger multi-pfaffian 
expansions than we have explored. As another test of 
the variational potential of pairing, we have employed 
the fully-antisymmetrized independent pairs wave func- 



FIG. 3: (color online) Same as Fig. [2] but for C dimer. 



tion in pfaffian form and we have found that it does not 
lead to additional gains in correlation energy. We there- 
fore conclude that more general functional forms together 
with more robust large-scale optimization methods might 
be necessary in order to obtain further improvements. 
The gains in correlation energy for pfaffians come from 
improved fermion nodes which are significantly closer to 
the exact ones than the HF nodes and exhibit the correct 
topology with the minimal number of two nodal cells. 

In the second part of the paper, we have presented 
the application of pfaffian and multi-determinantal wave 
functions with backflow correlations to chemical systems. 
Results for two testing cases of C atom and its dimer 
show promising gains in correlation energies, decreases in 
variances and improvements in the nodal structure. Our 
results also indicate that accurate description of molec- 
ular systems with fixed-node errors bellow 1% of cor- 
relation energy requires optimized multi-reference wave 
functions and inclusion of backflow then appears less fa- 
vorable considering significant computational cost, espe- 
cially for evaluation of the nonlocal pseudopotential op- 
erators. 
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APPENDIX A: PROOF OF CAYLEY'S IDENTITY given by 



In order to prove the statement in Eq. ([7]) we will pro- 
ceed by induction. For n = 2 it is true that 



det 



6ia 
-an 



= pf 



bia 
-bis 



pf 



a i2 
-ai2 



For even n greater than 2, determinant of the matrix of 
interest can be expanded through its cofactors as 





- 


b u 


bl3 ■ 


■ &i, n ' 






— Ol2 





023 ■ 


■ 0,2,n 




det 


-Ol3 


-023 


. 


■ «3 : n 


= ]T-ai, fe C(fc,l) 

k 




- — 0.1, n 


— 0.2,n 


— 0.3,1% ■ 


. o . 





^2^2-a 1<k bi,iC(k,l-AJ)- 



(Al) 



The cofactor can be written as 

C(k, 1; 1, I) = (-l) fe+ ' +1 det [A(k, 1; 1, I)} , (A2) 
where the cofactor matrix is given by 








a.23 


0.2, k 


■ 0.2,71 




— 0.23 





az,k 


■ ■ 03, n 


A(fc,i;i,0 = 


— ai,i 


— 0.3,! • 


■ —a,h,i ■ 


■ ■ ai, n 




—0,2,n 


— 0.3, n . 


■ — 0,k,n ■ 


■ 



(A3) 

At this point we would like to use an induction step and 
rewrite the determinant cofactor as a product of two pfaf- 
fians [Cayley's identity Eq. This would allow us to 
demonstrate that the expansion is identical to the expan- 
sion of pfaffians in minors. In order to do so, however, we 
have to shift the k-th column by pair column exchanges, 
so it becomes the last column and, similarly, we have to 
shift the l-th row by pair exchanges, so it becomes the 
last row. This involves k pair exchanges of columns and 
I pair exchanges or rows and can be represented by uni- 
tary matrices Uk and [//. It is necessary to invoke these 
operations so that the matrix gets into a form directly 
amenable for the Cayley's identity, i.e., the matrix has 
to be in a manifestly skew-symmetric form. (The sign 
change from the row/columns exchanges will prove irrel- 
evant as we will show below.) The transformed matrix is 
given by 



A\k,l;l,l) = U k A(k,l;l,l)Ut 



(A4) 



and has all zeros on the diagonal with the exception of 
the last element which is equal to — a k i- The last row is 



v r =(— a.2,1, • • ■ , — o,k-i,i, — ctfc+ijj, . . . 

. . . , — at-i t i, a/,j+i, ■ ■ • , ai t n, —au,i), (A5) 
while the last column is given as following 

T / 

v c — {a,2,k, ■ ■ ■ , Ok-i,kt —ak,k+i, ■ ■ ■ 

■ ■ ■ ) — flfc,J_l, — a,k,i+i, • • ■ , —Ok,rn ~ a k,i) T ■ (A6) 

The only non-zero diagonal clement —a k ,i can be elimi- 
nated, once we realize that its cofactor contains a deter- 
minant of a skew-symmetric matrix of odd degree which 
always vanishes (proof by Jacobi)^ 

Now we are ready to perform the induction step, 
namely to use the property that the determinant of a 
2(n— l)x2(n— 1) matrix can be written as given by the 
Cayley's identity Eq. ©. We obt am 

det[U k A(k, 1; l,l)U{\ = det[A'(fc, 1; 1, 1)] (A7) 
= P f[A'(l,A;;l,fe)]pf[A'(l,Z;l,0]. 

We can now apply the inverse unitary transformations 
and shift back the columns (and by the skew-symmetry 
the corresponding rows) in the first pfaffian and, simi- 
larly, the rows (and corresponding columns) in the sec- 
ond. This enables us to write 

P f[A'(i J *;i,*)]pf[^ , (i,i;i»0] 

= pf[U^A(l, k- 1, k)Ui] pf[U k A(l,l; 1, Ip- 1 ] 
= P f[A(l,fc;l,fc)]pf[A(M;l,OL (A8) 



where we have used the identity given by Eq. (|5d[) . We 
can therefore finally write 

C(k, 1; 1, 1) = (-l) fc+/+1 pf [A(l, k; 1, fc)]pf [A(l, I; 1,1)) 
= -Pc(ai,k)Pc(ai,i), (A9) 

where P c denotes a pfaffian cofactor as defined in Q. 
Therefore, the determinant expansion in Eq. (|A1[) equals 
to 

-ai t kh t iC(k, 1; 1,1) = ^ ai, k bi,iP c {ai, k )Pc(ai,i) 



k.l 



k.l 



pf[A]pf[B] 



(A10) 



with matrices A and B defined as in Eq. ([5]). This con- 
cludes the proof of the more general form of the Cayley's 
identity. Note, if B = A, we trivially obtain well-known 
formula for the square of pfaffian [Eq. (|5b[) ]. 



APPENDIX B 
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TABLE III: Total energies for C, N and O atoms and their dimers with amounts of the correlation energy recovered in VMC 
and DMC methods with wave functions as discussed in the text. Unless noted otherwise, the numbers in parentheses are the 
statistical errors in the last digit from corresponding QMC calculation. Energies are in Hartree atomic units. For C, N, O 
atoms we used the correlation energies by Dolg 55 (0.1031, 0.1303 , 0.1937 a.u.). For the estimation of correlation energies of 
dimers we needed accurate HF energies at experimental distances 5 - 6 - and the estimated exact total energies. Each exact total 
energy was estimated as a sum of total energies of constituent atoms minus experimental binding energ y 56,57 ' 58 adjusted for 
experimental zero-point energy.— 



Method 


WF 


C 


E corr [%] 


N 


E corr [%] 


o 


E corr [%] 


HF 


S 


-5.31471 





-9.62892 





-15.65851 





VMC 


SJ 


-5.3939(4) 


76.8(4) 


-9.7375(1) 


83.3(1) 


-15.8210(6) 


83.9(3) 




BCS 


-5.4061(2) 


88.6(2) 


-9.7427(3) 


87.3(2) 


-15.8250(3) 


86.0(2) 




STU 


-5.4068(2) 


89.3(2) 


-9.7433(1) 


87.8(1) 


-15.8255(3) 


86.2(2) 


DMC 


SJ 


-5.4061(3) 


88.6(2) 


-9.7496(2) 


92.6(2) 


-15.8421(2) 


94.8(1) 




BCS 


-5.4140(2) 


96.3(2) 


-9.7536(2) 


95.7(2) 


-15.8439(4) 


95.7(2) 




STU 


-5.4139(2) 


96.2(2) 


-9.7551(2) 


96.8(1) 


-15.8433(3) 


95.4(2) 


Est. 


Exact 


-5.417806 


100 


-9.759215 


100 


-15.85216 


100 


Method 


WF 


c 2 


E corr [%] 


N 2 


E CO rr[%] 


o 2 


E corr [%] 


HF 


S 


-10.6604 





-19.4504 





-31.3580 





VMC 


SJ 


-10.9579(4) 


72.9(1) 


-19.7958(5) 


80.0(1) 


-31.7858(6) 


79.6(1) 




BCS 


-11.0059(4) 


84.7(1) 


-19.8179(6) 


85.0(1) 


-31.8237(4) 


86.7(1) 




STU 


-11.0062(3) 


84.8(1) 


-19.821(1) 


85.8(2) 


-31.8234(4) 


86.6(1) 


DMC 


HF 


-11.0153(4) 


87.0(1) 


-19.8521(3) 


93.0(1) 


-31.8649(5) 


94.3(1) 




BCS 


-11.0416(3) 


93.5(1) 


-19.8605(6) 


94.9(1) 


-31.8664(5) 


94.6(1) 




STU 


-11.0421(5) 


93.6(1) 


-19.8607(4) 


95.0(1) 


-31.8654(5) 


94.4(1) 


Est. 


Exact c 


-11.068(5)° 


100.0(10) 


-19.8825(6)'' 


100.0(1) 


-31.8954(1)° 


100.0(1) 



"There is rather large discrepancy in the experimental values of 
C 2 binding energy (141.8(9)^ 143(3)S and 145.2(5) kcal/mol)S 
For the estimation of exact energy we have taken the average value 
of 143(3) kcal/mol. 

''Experimental binding energies taken from ref. |56|| . 

The error bars on estimated exact total energies are due to ex- 
periment. 
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TABLE IV: VMC and DMC energies and variances of local energy for Slater-Jastrow (SJ), pfaffian-Jastrow (PF) and CI- 
Jastrow (CI) trial wave functions with backflow (BF) correlations for C atom. Notation is the same as in Fig. [2] 



Method 


WF 








Ng 2 




E [a.u.] 


o 2 [a.u. 2 ] 


E corr [%] 


HF 


S 




_ 


_ 


_ 


_ 


-5.31471 


_ 


0.0 


VMC 


SJ 


- 


_ 


_ 


_ 


_ 


-5.3990(1) 


0.0677 


81.8(1) 




SJBF2B 


11 


22 


_ 


_ 


33 


-5.4011(2) 


0.0544 


83.8(2) 




SJBF3B 


- 


_ 


128 


128 


256 


-5.4023(3) 


0.0504 


85.0(3) 




SJBF23B 


4 


8 


128 


128 


268 


-5.4020(2) 


0.0498 


84.7(2) 




PF 


_ 


_ 


_ 


_ 


_ 


-5.4083(2) 


0.0626 


90.8(2) 




PFBF2B 


ii 


22 


_ 


_ 


33 


-5.4097(1) 


0.0427 


92.1(1) 




PFBF23B 


4 


8 


128 


128 


268 


-5.4107(1) 


0.0411 


93.1(1) 




cr 


_ 


_ 


_ 


_ 


_ 


-5.4127(1) 


0.0447 


95.0(1) 




CIBF2B 


11 


22 


_ 


_ 


33 


-5.4131(3) 


0.0427 


95.4(3) 




CIBF23B 


4 


8 


128 


128 


268 


-5.4140(1) 


0.0342 


96.3(1) 


DMC 


SJ 


_ 


_ 


_ 


_ 


_ 


-5.4065(3) 


_ 


89.0(3) 




SJBF2B 


11 


22 


- 


- 


33 


-5.4090(3) 


- 


91.5(3) 




SJB1 3h> 






128 


128 


256 


r ,1 AO C / 0\ 

-5.4085(3) 




91.0(3) 




SJBF23B 


4 


8 


128 


128 


268 


-5.4094(3) 




91.8(3) 




PF 












-5.4137(3) 




96.0(3) 




PFBF2B 


11 


22 






33 


-5.4145(3) 




96.8(3) 




PFBF23B 


4 


8 


128 


128 


268 


-5.4152(3) 




97.5(3) 




CI 












-5.4178(1) 




100.0(1) 




CIBF2B 


11 


22 






33 


-5.4177(3) 




99.9(3) 




CIBF23B 


4 


8 


128 


128 


268 


-5.4174(2) 




99.6(2) 


Est. 


Exact 












-5.417806 




100.0 



"Wave function consists of 100 determinants re-optimized in 
VMC. 



TABLE V: Slater-Jastrow (SJ), pfaffian-Jastrow (PF) and Cl-Jastrow (CI) wave functions with backflow (BF) correlations 
for C dimer. Notation is the same as in Fig. [21 



Method 


WF 








Ne 2 


N p 


E [a.u.] 


o 2 [a.u. 2 ] 


E corr [%] 


HF 


S 












-10.6604 




0.0 


VMC 


SJ a 












-10.9936(4) 


0.179 


81.7(1) 




SJBF2B 


11 


22 






33 


-11.0012(3) 


0.144 


83.5(1) 




SJBF23B 


4 


8 


128 


128 


268 


-11.0014(2) 


0.141 


83.6(1) 




PF b 












-11.0171(2) 


0.160 


87.4(1) 




PFBF2B 


11 


22 






33 


-11.0223(3) 


0.123 


88.7(1) 




PFBF23B 


4 


8 


128 


128 


268 


-11.0223(2) 


0.128 


88.7(1) 




cr 












-11.0420(4) 


0.112 


93.6(1) 




CIBF2B 


11 


22 






33 


-11.0440(3) 


0.100 


94.0(1) 




CIBF23B 


4 


8 


128 


128 


268 


-11.0438(3) 


0.123 


94.0(1) 


DMC 


SJ 












-11.0227(2) 




88.8(1) 




SJBF2B 


11 


22 






33 


-11.0269(4) 




89.9(1) 




SJBF23B 


4 


8 


128 


128 


268 


-11.0280(3) 




90.1(1) 




PF 












-11.0419(9) 




93.5(2) 




PFBF2B 


11 


22 






33 


-11.0443(6) 




94.1(2) 




PFBF23B 


4 


8 


128 


128 


268 


-11.0447(3) 




94.2(1) 




CI 












-11.0579(5) 




97.5(1) 




CIBF2B 


11 


22 






33 


-11.0580(4) 




97.5(1) 




CIBF23B 


4 


8 


128 


128 


268 


-11.0585(5) 




97.7(1) 


Est. 


Exact 












-11.068(5) 




100.0 



a Slatcr determinant contains PBE DFT orbitals. 
b Same PBE DFT orbitals are used also in PF wave function. 
c Uses natural orbitals with weights of the 148 determinants re- 
optimized in VMC. 
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